What’s New in Probability & Statistics

Devendra Kapadia 
Wolfram Research, Inc.

Outline of the Talk

Time Series Operations and Analysis

Transformations of Random Processes

Hidden Markov Models

Probability and Statistics in Mathematica 10

Time Series Operations and Analysis

TimeSeries and EventSeries

Mathematica 10 introduces TimeSeries and EventSeries for representing continuous and discrete time series, respectively.

For example, the following is the temperature variation for Houston, Texas on October 16:

ts=TimeSeries[];

DateListPlot[ts, PlotTheme -> "OpenMarkers"]

List the temperatures at different times during the day:

ts["Values"]

Find the temperature at a specific time:

ts["Oct 16, 13:17:35"]

Personal Fitness Analysis

This time series contains the number of steps taken daily by a person during a period of five months.

StepData=TimeSeries[];

Plot the daily step counts together with a healthy goal of 10,000 steps per day.

tenKlevel = {{{2013, 04, 01}, 10^4}, {{2013, 09, 01}, 10^4}};
opt = {Joined -> {False, True}, Filling -> 0, ImageSize -> Medium}; DateListPlot[{StepData, tenKlevel}, opt]

Compute the basic statistics for the daily number of steps.

Grid[Transpose@Map[{#, Floor[#[StepData]]} &, {Min, Mean, Median, Max}], Alignment -> {Left, Center}, Dividers -> All, Spacings -> {1.5, 1.}]

Analyze the frequency of daily counts.

Histogram[StepData, 10, ImageSize -> 300]

Operations on Time Series

A wide variety of functions are available for operating on time series e.g. TimeSeriesWindow, TimeSeriesThread, etc.

For example, we can resample an irregular time series using TimeSeriesResample as follows:

irrts = TimeSeries[RandomFunction[PoissonProcess[0.73], {0, 30}]]

regts = TimeSeriesResample[irrts, 0.5, ResamplingMethod -> {"Interpolation", InterpolationOrder -> 0}]

ListPlot[{regts, irrts}, Filling -> {2 -> Axis}, Joined -> {False, False}, FillingStyle -> Thick, PlotTheme -> {"Scientific", "PlotMarkers"}, PlotLegends -> {"resampled time series", "original time series"}, ImageSize -> 300, Ticks -> False]

Time Series Analysis

Number of mink furs sold yearly by the Hudson Bay Company, Canada from 1850 to 1911:

sales=TimeSeries[];

ListPlot[sales, Joined -> True, PlotRange -> All]

Fit a model to the data:

saleModel = TimeSeriesModelFit[sales]

Simulate the values of the model for the next 10 years:

sim = RandomFunction[saleModel, {0, 10}, 150]

ListPlot[sim, Joined -> True, PlotRange -> All]

The simulated values may be negative, but the mean function provides a positive forecast:

forecast = TimeSeriesThread[Mean, sim];

Min[forecast]

Plot the data with the forecast given by the simulation:

Show[ListLinePlot[sales], ListLinePlot[sim, PlotStyle -> Directive[Opacity[.2]]],
 ListLinePlot[forecast, PlotStyle -> Black], PlotRange -> All]

Transformations of Random Processes

TransformedProcess

TransformedProcess allows the user to define new processes as point transformations of existing random processes:

TransformedProcess[expr, x\[Distributed] proc,t]

It can be used to gain insight into stochastic calculus, finance, insurance, and other fields.

The Stochastic Exponential Function

Define the stochastic exponential of the Wiener process as a transformed Wiener process.

ℰw=TransformedProcess[(E^(w[t]-(((σ^2) t)/2))),w\[Distributed]WienerProcess[0,σ],t];

Simulate the process.

data=RandomFunction[ℰw/.σ->1,{0,10,0.01},4];

ListLinePlot[data,PlotRange->All]

Compute the mean and variance for a time slice of the process.

{Mean[ℰw[t]],Variance[ℰw[t]]}

Compare with the result using the corresponding SDE du(t)==u(t) db(t).

ℰsde=ItoProcess[\[DifferentialD]u[t]==u[t]\[DifferentialD]b[t],u[t],{u,1},t,b\[Distributed]WienerProcess[0,σ]];

{Mean[ℰsde[t]],Variance[ℰsde[t]]}

Define the stochastic exponential of the compensated Poisson process as a transformed Poisson process.

ℰp = TransformedProcess[2^p[t] Exp[-λ t], p \[Distributed] PoissonProcess[λ], t];

ListLinePlot[RandomFunction[ℰp/.λ->1,{0,10,0.01},4],PlotRange->All]

{Mean[ℰp[t]],Variance[ℰp[t]]}

A Jump-Diffusion Model

Define a Merton-type jump-diffusion model for option pricing.

Merton𝒫=TransformedProcess[(v-((σ^2)/2)) t+σ w[t]+j[t],{w\[Distributed]WienerProcess[r,σ],j\[Distributed]CompoundPoissonProcess[λ,NormalDistribution[μ,δ]]},t];

Simulate the process.

data=RandomFunction[Merton𝒫/.{v->1.2,σ-> 0.7,λ->1.3,μ->0.92,δ->0.425,r->1},{0,3,0.02},3];

ListLinePlot[data]

Compute slice properties for the process.

{Mean[Merton𝒫[t]],Skewness[Merton𝒫[t]],Kurtosis[Merton𝒫[t]]}

Approximate the slice distribution from a sample.

sliceSample=RandomVariate[Merton𝒫[1]/.{v->1.2,σ-> 0.7,λ->1.3,μ->0.92,δ->0.425,r->1},10^6];

SmoothHistogram[sliceSample]

Hidden Markov Models

HiddenMarkovProcess

A hidden Markov model (HMM) is a statistical Markov model in which the system being modeled is assumed to be a Markov process with hidden states.

HiddenMarkovProcess provides functionality for simulation, estimation of parameters, finding the hidden states in such a model, support for silent states, etc.

An Unfair Casino

A dealer at a casino uses a fair coin or a biased coin (75% of heads), and switches the coin used 2 out of 10 times.

The unfair casino may be represented using HiddenMarkovProcess as follows:

unfairCasino=HiddenMarkovProcess[{(1/2),(1/2)},{{(8/10),(2/10)},{(2/10),(8/10)}},{{(1/2),(1/2)},{(3/4),(1/4)}}];

For comparison, we also define the corresponding fair casino.

fairCasino=HiddenMarkovProcess[{1},{{1}},{{(1/2),(1/2)}}];

The following plot compares the distribution of the head counts for the unfair and fair casinos.

HeadsCount[td_TemporalData]:=Total[2-td["ValueList"],{2}]

Histogram[{HeadsCount[RandomFunction[unfairCasino,{0,29},3000]],
HeadsCount[RandomFunction[fairCasino,{0,29},3000]]},Automatic,"Probability",
ChartLegends->{"Unfair casino","Fair casino"},
PlotTheme->"Detailed",
PlotLabel->"Distribution of head counts in a 30 coin throws"]

We now generate a sequence of observed tossings for the unfair casino.

alphabet=Association["head"->1,"tail"->2];

{PDF[unfairCasino[t],alphabet["head"]],PDF[unfairCasino[t],alphabet["tail"]]}

Short[tosses=RandomChoice[{5/8,3/8}->{"head","tail"},300]]

With the assumed model of dealer behavior, we can find the most probable sequence of underlying hidden states.

Short[FindHiddenMarkovStates[alphabet/@tosses,unfairCasino,"ViterbiDecoding"]]

ListPlot[FindHiddenMarkovStates[alphabet/@tosses,unfairCasino,"ViterbiDecoding"],ImageSize->Medium,Filling->Axis,Ticks->{Automatic,{{1,"Fair"},{2,"Unfair"}}},PlotLabel->"Viterbi decoding"]

Probability and Statistics in Mathematica 10

A Brief History

The following is a brief history of Probability and Statistics development over the last three versions of Mathematica.

Version 8

Random variate framework introduced [solvers] [estimation & testing] [distributions] [derived and non-parametric]

Version 9

Random processes framework defined [slices, parametric processes, estimation, time series, SDE ] [Markov chains and queues ]

RLink` introduced [web]

Random variate framework refined and extended [overview] [reliability analysis] [survival analysis]

Version 10

TimeSeries + EventSeries data containers and their support in the system [time-series]

Automation + enhancements in time series modeling [time-series-modeling]

New processes + better integration with P&S framework [enhanced-random-processes]

HiddenMarkovProcess + FindHiddenMarkovStates

TransformedProcess

GARCHProcess + ARCHProcess

Comprehensive support for computation with slices

Better intergration of processes in Probability, Expectation, etc. solvers.

RLink` can work with arbitrary custom installation of R.

Documentation

There are now around 30,000 examples for Probability and Statistics alone.

perm = Reverse[{1, 3, -1, 4, 5, 2, -3, -2}];
data = {{1820, 4026, 6255, 1143, 3623, 470, 158, 113}, {2482, 5402, 9628, 1693, 4735, 614, 188, 125}, {2641, 5987, 11727, 1863, 5743, 777, 632, 124}}; legendLabels = {"Basic Examples", "Properties & Relations", "Scope", "Options", "Applications", "Possible Issues", "Neat Examples", "Generalizations & Extensions"};
docsBC = BarChart[Part[data, All, perm], ChartLayout -> "Stacked", ChartLegends -> Part[legendLabels, perm], ChartLabels -> {{"v8", "v9", "v10"}, None}, ChartStyle -> 98, PlotLabel -> "Section Counts", ImageSize -> 300];
CellPrint[ExpressionCell[docsBC, "Output", Background -> None, ShowCellLabel -> False, CellFrame -> False]]

The documentation was greatly improved in Mathematica 10 to address the needs of beginners as well as common issues raised by our users, and to provide more neat examples for the functionality.

Summary

Mathematica 10 offers a comprehensive system for doing symbolic/numerical computation, data analysis, visualization and programming for a wide variety of tasks in Statistics.

The Statistics Team

Thank you from us all!